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ABSTRACT 


Previous  studies  have  shown  that  terminal  velocity  raindrops  striking  a  smooth  water 
surface  create  oscillating  bubbles  that  radiate  significant  underwater  sound  energy.  Those 
studies  identified  two  diameter  ranges  that  produce  bubbles:  small  drops  (0.8-1. 1  mm 
diameter)  which  produce  bubbles  by  one  mechanism  and  large  drops  (2.2-4.6  mm 
diameter)  which  create  bubbles  by  a  different  mechanism.  Effects  of  oblique  incidence 
have  been  studied  only  for  small  drops.  Average  energy  spectra  were  calculated  for  a 
range  of  raindrop  sizes  striking  a  smooth  water  surface. 

This  work  deals  with  the  real  life  situation  of  large  raindrops  of  a  size  often  present  in 
heavy  rainfall  (4.6  mm  diameter)  striking  a  sloped  water  surface.  Terminal  velocity  is  used 
to  simulate  natural  rainfall,  and  the  sloped  surface  is  used  to  simulate  the  surface  gravity 
waves  of  a  natural  sea.  The  effects  of  a  sloped  water  surface  on  the  frequency  spectra  and 
energy  for  4.6  mm  raindrops  are  estimated. 

By  comparing  energy  spectra  generated  by  single  drops  in  an  anechoic  laboratory  tank 
to  underwater  sound  spectra  measured  at  sea,  it  will  be  possible  to  estimate  heavy  rainfall 
rate  by  means  of  remote  underwater  listening  devices. 
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I.  INTRODUCTION 


One  of  the  goals  of  the  Naval  Postgraduate  School  (NPS)  Raindrop  Lab  is  to 
understand  the  physics  of  underwater  sound  production  by  natural  rainfall.  Previous 
studies  have  shown  that  significant  underwater  sound  is  generated  by  the  impact  of 
raindrops  upon  a  water  surface  and  even  more  so  by  oscillating  bubbles,  when  formed 
[Franz,  1959],  Once  the  sound  pressure  is  measured,  the  energy  spectrum  of  a  single  drop 
can  be  calculated.  Figure  1 . 1  shows  a  typical  time  signal  of  an  impact  of  a  large  diameter 
(4.6  mm )  drop  which  produced  two  bubbles  .  The  corresponding  energy  spectrum  is  also 
shown. 

By  cataloging  the  sound  produced  by  each  raindrop  size,  it  is  possible  to  predict  the 
sound  produced  by  natural  rainfall  [Nystuen  et  al  1992].  This,  in  turn,  will  allow  rainfall 
to  be  monitored  using  passive  listening  devices. 

Raindrops  can  be  categorized  by  size  and  acoustic  signal  [Medwin  et  al,  1992].  Table 
1.1  summarizes  the  common  drop  size  categories  and  their  corresponding  sources  of 
sound. 


TABLE  1.1.  Raindrop  Sizes  and  their  Sources  of  Sound. 


Nomenclature 

Rundrop  Diameter  (mm) 

Sources  of  Underwater  Sound 

Minuscule 

0-0.8 

Impact 

Small  (Type  I) 

0.8-1. 1 

Impact  and  Bubbles 

Mid'Size 

1. 1-2.2 

Impact 

Large  (Type  II) 

2.2  and  larger 

Impact  and  Bubbles 

1 


(Joules  /  H?] 


freg  [kHz] 

Figure  1.1.  Typical  graph  of  a  time  series  of  on-a.\is  acoustic  pressure  at  1  m  for  an  impact  and  two 
resulting  bugles  (upper  plot)  and  the  corresponding  energy  spectrum  (lower  plot). 
{Sampling  frequency  50  kHz.;  4.6  mm  drop}. 
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The  most  important  size  categories  are  small  and  large  drops  as  these  two  sizes 
consistently  produce  oscillating  bubbles  which  are  the  dominant  sound  source.  (The  Type 
I  and  Type  II  bubble  formation  mechanisms  are  described  in  Chapter  III.) 

All  experiments  at  NPS  have  been  conducted  with  raindrops  at  terminal  velocity. 
Previously,  however;  the  large  drops  were  impacting  at  normal  incidence  upon  smooth 
water  surfaces.  A  normal  angle  of  incidence  does  not  describe  a  natural  sea  state. 
Although  these  experiments  were  outstanding  for  determining  the  basic  physics  of 
underwater  sound  production,  a  more  realistic  scenario  is  required. 

Kurgan  [1989]  conducted  experiments  for  small  drops  at  terminal  velocity  and  at 
various  angles  of  incidence.  The  influence  of  angle  of  incidence  on  small  drops  was  very 
significant.  His  non-normal  incidence  angles  were  generated  by  placing  a  fan  above  the 
point  of  impact  and  forcing  the  small  drops  to  strike  the  smooth  surface  at  an  oblique 
angle.  Logistically  large  raindrops  can  not  be  significantly  deflected  using  a  fan. 

The  purpose  of  this  thesis  is  to  apply  the  theory  and  knowledge  gained  from  past 
experiments  with  smooth  surfaces  to  large  raindrops  impacting  upon  a  sloped  water 
surface.  This  simulates  ocean  surface  gravity  waves.  Large  drops,  with  a  diameter  of  4.6 
mm  at  terminal  velocity,  were  used  throughout  the  experiment.  The  results  will  be 
compared  to  those  for  smooth  surface  experiments. 


3 


II.  LABORATORY  FACILITIES 


A.  EXPERIMENTAL  SETUP 

1.  Drop  Shaft 

The  laboratory  facilities  available  for  this  experiment  are  unique  to  the  Naval 
Postgraduate  School.  A  utility  shaft  with  a  height  of  26  meters  (and  cross  section  of  3  m 
X  3  m)  empties  into  a  room  containing  a  1.5  m  diameter,  redwood  lined  anechoic 
cylindrical  tank  .  [Figure  2.1]  Installed  in  the  shaft  is  a  9"  diameter  plastic  tube  that 
reduces  the  amount  of  interference  from  air  drafts.  The  26  m  height  allows  drops  of  all 
sizes  normally  found  in  natural  rainfall  to  reach  terminal  velocity. 
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2.  Intravenous  (IV)  Dropper 

A  standard  medical  IV  bag  with  surgical  tubing  is  used  to  generate  large  raindrops. 
Attached  to  the  tubing  is  a  tip  calibrated  to  produce  50  p.L  drops  (4.6  mm  diameter).  The 
accuracy  of  the  tip  was  verified  by  measuring  the  volume  of  100  drops  six  times.  Each 
repetition  was  ±  5  %  of  the  volume  expected  for  100  drops. 

3.  Anti-Wind  Tube 

New  to  the  lab  is  the  anti-wind  tube.  It  is  constructed  of  thirty-three  needlepoint 
hoops  (  nine  inch  diameter)  and  two  50  ft.  lengths  of  typical  plastic  liner  for  gardens.  The 
two  lengths  are  overlapped  by  13  feet.  The  hoops  are  equally  spaced  the  entire  length  of 
the  tube.  The  inner  ring  of  the  hoop  is  rolled  inside  of  the  plastic.  The  outer  ring  is  made 
tight  around  the  inner  ring.  This  forms  the  plastic  into  a  9"  x  87'  cylinder. 

The  addition  of  the  tube  increases  the  accuracy  of  the  impact  of  the  raindrop  in 
relation  to  the  location  of  the  hydrophone.  Previous  work  measured  impacts  with 
estimates  of  the  random  horizontal  distances  from  the  hydrophone  (within  20  cm). 
Presently,  we  are  consistently  able  to  impact  the  surface  of  the  water  within  4  cm  of  the 
epicenter  of  the  underwater  hydrophone.  The  increased  accuracy  reduces  the  amount  of 
range  and  angle  correction  required. 

4.  Anechoic  Tank 

In  the  room  where  the  anti-wind  tube  ends  is  a  cylindrical  redwood  tank  with  a 
diameter  and  a  height  of  1.5  m.  The  tank  is  made  anechoic  with  a  lining  of  redwood 
wedges.  The  tank  contains  filtered  salt  water  acquired  from  the  Monterey  Bay  Aquarium. 
Because  data  sets  were  taken  on  different  days,  the  salinity  was  measured  for  each  data 
set. 

The  tank  houses  the  wave  generator  fi-ame  and  an  LC-10  hydrophone. 


5 


5.  Wave  Generator  Frame  and  Motor 

Also  new  to  the  lab  is  the  wave  generator  frame  and  motor.  It  is  inserted  into 
the  tank  and  resembles  a  box  kite.  It  is  constructed  of  a  redwood  frame  that  supports 
walls  made  of  Mylar.  Three  sides  are  stationary  while  the  fourth  Mylar  wall  is  hinged  at 
the  bottom  to  allow  for  a  paddle  motion.  The  paddle  is  attached  to  a  motor  mounted  level 
with  the  top  of  the  tank.  The  dimensions  of  the  frame  are  26"  x  26"  x  4 1"  (L  x  W  x  H). 

The  wave  generator  remains  in  the  tank  for  experiments  involving  both  smooth  and 
rough  surfaces.  Mylar  was  used  because  of  its  anechoic  properties.  It  allows  sound  to 
pass  but  reflects  surface  water  waves.  This  maintains  the  anechoic  integrity  of  the  tank 
and  still  allows  for  the  creation  of  a  roughened  surface. 

The  motor's  shaft  generates  circular  motion.  This  motion  is  translated  to  "linear" 
motion  to  drive  the  paddle  via  a  mechanical  eccentric  coupling. 

6.  Video  Camera 

A  Sony  8  mm  video  camera  placed  level  with  the  water  surface  is  used  to  film  the 
sloped  surface  experiment.  Extracted  from  the  videotape  are  the  slope  of  the  surface  and 
the  depth  of  the  hydrophone  at  the  time  of  impact.  The  camera  speed  is  30  frames  per 
second. 

To  assist  with  the  determination  of  the  slope,  a  grid  of  vertical  strings  is  suspended 
in  the  tank  .  The  grid  is  filmed  in  the  background  while  the  drops  impact  the  surface.  The 
spacing  between  vertical  lines  is  2  cm.  A  reference  line  is  marked  on  the  grid  before  the 
surface  is  roughened.  This  allows  for  the  measurement  of  the  depth  of  the  hydrophone. 
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7.  Hydrophone 

An  LC-10  hydrophone  [Figure  2.2]  from  Celesco  Transducer  Products,  Inc.  is 
used  to  measure  the  sound  pressure  of  the  impacts  and  oscillating  bubbles.  It  is  positioned 
at  a  depth  of  6  cm  for  both  types  of  experiments.  For  the  sloped  surface  experiment  a 
video  camera  is  used  to  determine  the  instantaneous  depth  of  the  LC-10.  The  hydrophone 
is  suspended  by  three  supports  of  fishing  line.  The  supports  are  positioned  at  intervals  of 
120°  to  minimize  the  amount  of  movement  caused  by  a  passing  wave. 


s!a 

ibbbEI^I 

Figure  2.2.  LC-10  Hydrophone  (All  dimensions  are  in  inches). 


8.  Amplifiers  and  Filters 

The  signal  from  the  hydrophone  is  connected  to  an  Ithaco  1201  Pre-amp  set  at  a 
gain  of  100.  The  signal  is  then  connected  to  a  Krohn-Hite  3202R  band  pass  filter.  The 
frequencies  passed  are  between  1  kHz  and  30  kHz  in  both  pieces  of  equipment. 

9.  ComputerScope® 

The  signal  is  then  connected  to  an  IBM  clone  286  computer  with  an  analog  to 
digital  converter  (A/D).  The  digital  data  acquisition  card  and  its  software, 
ComputerScope  (sold  by  RC  Electronics),  is  used  for  all  data  acquisition. 

Amplitude  resolution  was  twelve  bits.  The  sampling  frequency  used  was  50  kHz. 
This  allowed  for  a  time  series  of  320  ms  duration.  Previous  work  was  sampled  at  125 
kHz  and  250  kHz.  Their  respective  time  series  lengths  were  128  ms  and  64  ms. 
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III.  BACKGROUND 


A.  BUBBLE  FORMATION  MECHANISMS 

Two  distinct  bubble  formation  mechanisms  have  been  identified  [Snyder,  1990],  The 
Type  I  mechanism  pertains  to  small  diameter  raindrops  (0.8  mm  to  1.1  mm)  whereas  the 
Type  II  mechanism  pertains  to  large  diameter  raindrops  (>  2.2  mm).  At  terminal  velocity 
the  raindrop  diameter  ranges  of  d  <  1.1  mm  and  1.1  mm  <  d  <  2.2  mm  do  not  produce 
bubbles. 

1.  Type  I 

Small  raindrops  at  terminal  velocity  produce  Type  I  bubbles  100  %  of  the  time 
when  striking  a  smooth  surface  at  a  normal  angle  of  incidence.  As  the  angle  of  incidence 
increases  to  20®,  the  bubble  formation  percentage  drops  to  10  %  [Kurgan,  1989].  A 
Type  I  bubble  is  formed  when  the  base  of  a  conical  splash  crater  is  pinched  oflf  [Longuet- 
Higgins,  1990].  The  resulting  bubble  resonates  at  approximately  15  kHz.  Additional 
references  for  the  Type  I  mechanism  are  Pumphrey  et  al.^  [1989]  and  Oguz  and 
Prosperetti,  [1990],  This  mechanism  can  be  used  to  explain  the  sound  produced 
underwater  by  light  rain  [Pumphrey  et  al,  1989]  and  the  influence  of  wind  upon  the  sound 
produced  underwater  by  light  rain  [Nystuen,  1977]. 

2,  Type  II 

The  Type  II  mechanism  is  important  when  large  raindrops  are  present  in  the  rain. 
Large  raindrops  are  prevalent  during  heavy  rain. 

Type  n  bubbles  are  not  formed  as  consistently  as  Type  I  bubbles.  At  terminal 
velocity  and  normal  incidence  upon  a  smooth  water  surface,  it  has  been  shown  that  4.6 
mm  drops  form  at  least  one  bubble  50-65  %  of  the  time.  When  Type  n  bubbles  are 
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formed  they  resonate  between  1.6  and  10  kHz  [Jacobus,  1991;  Ostwald,  1992]  depending 
on  drop  diameter.  For  4.6  mm  drops,  the  dominant  resonance  frequencies  are  between  1.6 
and  2.0  kHz. 

Snyder,  Jacobus,  and  Ostwald  describe  Type  II  bubbles.  They  define  a  primary 
bubble  as  the  bubble  that  produces  the  largest  peak  to  peak  voltage  (pressure)  for  that 
drop  when  produced  by  the  Type  II  mechanism  (usually  60-70  ms  after  impact). 
Secondary  bubbles  were  defined  as  anything  else.  We  have  found  that  this  may  not  be 
appropriate.  It  would  be  better  to  define  the  "primary  bubble"  as  the  resulting  bubble  that 
contains  the  most  energy.  This  definition  is  independent  of  the  bubble  production 
mechanism.  Recent  experiments  with  large  drops  have  shown  that  some  previously 
defined  "primary"  bubbles  may  have  been  produced  by  the  Type  I  mechanism  generated  by 
late-arriving  aerosols  resulting  from  the  impact.  The  long  time  after  impact  (>100  ms) 
when  a  newly  defined  primary  bubble  begins  to  resonate  leads  us  to  believe  that  it  is  most 
likely  generated  by  a  Type  I  mechanism  (generated  by  a  small  drop,  an  aerosol,  striking 
the  surface).  Chapter  IV  addresses  the  significant  effect  of  a  sloped  surface  on  these  "late" 
bubbles. 

Figure  3.1  (frames  1-10)  shows  the  sequence  of  events  that  lead  to  the  formation 
of  a  Type  II  bubble.  This  phenomenon  was  first  observed  by  Snyder  [1990].  Slow 
motion  and  stop  motion  photography  of  400  frames  per  second  allowed  Snyder  the 
opportunity  to  sketch  the  sequence  . 

Frames  1  and  2:  A  flattened  raindrop  at  terminal  velocity  impacts  a  smooth 

water  surface  and  begins  to  form  a  canopy.  The  splash  generates 
numerous  aerosols. 

Frames  3  and  4:  The  canopy  continues  to  form  above  the  splash  crater. 

Frames  5, 6,  and  7;  Water  continues  to  flow  up  the  sides  of  the  canopy  and  the 

convergence  of  water  forms  upward  and  downward  moving 
turbulent  jets. 
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Frames  8  and  9:  The  downward  moving  jet  plunges  through  the  bottom  of 

the  crater. 

Frame  10:  This  jet  contains  entrained  air.  A  buoyant  force  causes  the 

entrained  air  to  break  off  and  form  a  resonating  bubble. 
This  action  occurs  only  if  the  downward  moving  jet  is 
canted. 
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B.  ENERGY  EQUATIONS 

1.  Conversion  of  Hydrophone  Voltage  to  Pressure 

To  calculate  energy,  the  hydrophone  voltage  signal  must  first  be  converted  to  a 
pressure  signal.  This  is  done  using  Equation  3.1, 


where  ^hyd  **  the  pressure  at  the  hydrophone  [Pa],  v  is  the  hydrophone  voltage  [volts], 
G  is  the  amplifier  gain,  and  ML  is  the  hydrophone  sensitivity  level  [volts/Pa]. 


2.  Corrections  to  Pressure  Signal 

a.  Correction  to  One  meter  on  Axis 

i.  Dipole  Radiation  Pattern.  Kurgan  [1989]  showed  that  the  bubble 
oscillation  pattern  was  that  of  a  dipole  (cos  0  term).  The  geometry  of  the  problem  is 
shown  in  Figure  3.2. 


ii.  Spherical  Spreading.  A  correction  is  applied  to  account  for  the 


1/r  divergence. 

Therefore,  all  pressures  measured  at  range  r  and  depth  z  were  converted  to 
1  meter  on  axis  using  Equation  3.2  where  to  =  1  m  and  cos  0  =  z/r. 


Pi 


m  axis 


Phyd  r_ 

COS  0  Tn 


(3.2) 
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6.  Near  Field  Correction  (NFC) 

Because  the  hydrophone  was  within  a  fraction  of  a  wavelength,  an  error 
between  the  measured  pressure  field  and  the  pressure  corrected  to  1  meter  on  axis  can 
occur.  This  can  be  corrected  in  the  frequency  domain  using  Equation  3.3  [Med^^  and 
Beaky,  1989], 


NFC  = 


(3.3) 


where  k  is  the  wave  number  (2idX  [m  ‘^])  and  r  is  the  range  [m]  as  shown  in  Figure  3.2. 

c.  Percentage  Correction  Factor 

The  voltage  time  series  of  a  single  drop  was  recorded  if  it  produced  at  least 
one  bubble.  Then  the  energy  was  calculated  and  averaged  for  all  drops  in  the  data  set. 
However,  because  only  data  from  drops  producing  bubbles  were  taken  and  the  fact  that 
raindrops  do  not  produce  Type  11  bubbles  100  %  of  the  time,  the  average  spectral  density 
was  obttuned  by  multiplying  the  calculated  energy  spectrum  by  a  percentage  correction 
factor  equal  to  the  percentage  of  drops  which  produced  bubbles  (determined  separately). 
Jacobus  and  Ostwald  extracted  the  impact  and  resonating  bubbles  from  each  time  series. 
They  then  calculated  the  average  energies  due  to  impacts  only,  and  then  repeated  the 
process  for  bubbles  only.  Jacobus  showed  that  the  impact  energies  are  significantly  less 
than  bubble  energies  (approximately  four  orders  of  magnitude).  To  obtain  the  "average" 
spectral  energy  due  to  bubbles,  the  percentage  correction  factor  was  applied  to  the 
calculated  bubble  spectrum. 

This  work  does  not  separate  the  impacts  and  bubbles  fi'om  the  time  series 
and  applies  the  percentage  correction  factor  to  the  energies  calculated  from  the  entire 
series.  The  error  in  this  approach  is  considered  insignificant. 
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d.  Summary 

Equation  3.4  shows  all  correction  factors  (except  the  percentage  correction 
factor)  considered  at  the  same  time  for  the  pressure  series. 


Plmmcu^fa^U  ^  |v(/)|  (3-4) 


V(f)  is  the  Fourier  transform  of  the  voltage-time  series.  Equation  3.4  combines  equations 
31,3.2,  and  3.3.  The  percentage  correction  factor  is  later  applied  to  the  ensemble  energy 
density,  rather  than  to  each  individual  drop  saved  as  data. 

3.  Spectral  Analysis  and  Total  Energy 

A  quantity  proportional  to  the  energy  density  spectrum  can  be  calculated  using 
Equation  3.5  [Ostwald,  1992]. 


£(/)  = 


Ndf 


'I  m  on  axis  farfteld ' 


Pa^  s 


J>0 


(3.5) 


N  is  the  number  of  points  in  the  FFT  and  the  frequency  resolution  is  given  by  df=l/(N  At). 

To  convert  the  units  to  Joules/Hz,  Equation  3.6  was  used  [Ostwald,  1992].  This 
is  the  energy  spectral  level  as  a  function  of  frequency,  and  a  subscript  is  added  to  avoid 


confusion  with  Equation  3.5. 


'Ite 

Ef=^  E(J) 


(3.6) 


This  is  the  evaluation  of  the  energy  density  as  shown  in  the  bottom  half  of  Figure  1.1. 


13 


The  total  energy  is  calculated  using  Equation  3.7.  The  subscript  "s"  is  used  to 
designate  a  spectral  calculation  and  to  avoid  confusion  with  the  energy  (£5^  which  is 
calculated  in  the  next  section  by  using  the  measured  theoretical  damping  constant  (5)  and 
the  peak  pressure  in  the  time  domain. 


/max 

I  Efdf 
/=0 


(3.7) 


4.  Energy  Calculation  using  the  Theoretical  Damping 
Constant  and  Peak  Pressure 

A  combination  of  equations  derived  by  Kurgan  [1989]  and  Scofield  [1992]  offers 
another  method  to  calculate  the  energy  of  an  individual  bubble  at  a  given  frequency  from 
the  peak  axial  pressure  and  the  damping  constant.  The  combination  is  derived  assuming 
dipole  radiation  and  the  theoretical  damping  constant  for  a  given  bubble  size.  Previously, 
Kurgan  empirically  calcualted  the  energy  and  did  not  consider  the  damping  constant. 
Scofield  was  interested  only  in  the  damping  constant  and  did  not  consider  the  dipole 
nature  of  a  bubble.  From  Figure  6.3.1  of  Clay  and  Medwin  [1977],  one  can  show  that  for 
the  frequencies  of  interest  for  this  work,  the  damping  constant  (5)  can  be  approximated  by 
5  =  .0025  f  (1^)  where  the  frequency  is  in  Hz.  The  energy  density  on  axis  from  an 
acoustic  dipole  [Scofield,  1992]  can  be  written  as 


Pi  m  on  axis 

S 


(3.8) 


where  the  pressure  is  the  PEAK  pressure  of  the  bubble  time  series. 
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The  energy  can  be  estimated  by  integrating  Equation  3.8  (using  a  dipole  radiation  pattern) 
[Kurgan,  1989], 

Eg  =  jj  Ea  cos^  0  dA  (3.9) 

The  elemental  area  is  dA  =  (2jc  r  sin  0  )(r  d0).  Evaluation  of  Equation  3.9  yields  the 
energy  equation  using  the  measured  theoretical  damping  constant  and  peak  pressure; 


E.= 


_  Pi  m  on  axis 


6poCfo  S 


[J] 


(3.10) 


IV.  RESULTS 


Smooth  surface  data  was  taken  as  well  as  sloped  (roughened)  surface  data.  It  was 
observed  that  a  sloped  surface  has  a  slightly  higher  percentage  of  at  least  one  bubble  being 
produced  (70-73  %)  than  a  smooth  surface  (50-65  %)  for  4.6  mm  diameter  drops.  The 
roughened  surface  data  was  analyzed  as  one  data  set  and  then  again  as  three  separate  data 
sets  broken  into  slope  categories.  A  data  set  consists  of  a  number  of  independent  drops, 
falling  at  terminal  velocity,  and  striking  either  a  smooth  or  a  sloped  surface.  A  drop  time 
series  was  saved  for  analysis  if  it  produced  at  least  one  bubble.  Because  our  time  series 
extended  further  than  ever  before  (320  ms),  many  more  bubbles  occurring  later  in  time 
were  detected.  Jacobus'  work  only  extended  to  64  ms  and  Ostwald's  time  series  ended  at 
128  ms.  The  duration  of  the  series  was  governed  by  the  available  computer  memory  and 
the  sampling  rate  of  the  data  acquisition  software  (ComputerScope®).  The  sampling  rates 
for  Jacobus  and  Ostwald,  respectively,  were  250  kHz  and  125  kHz.  This  work  uses  a 
sampling  rate  of  50  kHz.  This  rate  is  fast  enough  to  avoid  ambiguity  in  the  signals  with 
frequencies  of  interest  and  allows  for  a  longer  time  series.  The  following  presentation  of 
our  data  will  clearly  show  what  we  call  the  effects  of  a  sloped  water  surface  on  the 
underwater  sound  radiation  caused  by  large  raindrops.  For  this  thesis  only  the  largest 
drop  size  most  likely  to  occur  in  a  natural,  heavy  rainstorm  was  used  (4.6  mm  diameter). 

All  energies  were  calculated  using  MATLAB,  a  matrix  manipulation  program  by 
The  MathWorks,  Inc.  Appendix  A  contains  the  programs  used  to  gather,  process,  and 
plot  all  data.  The  spectral  analysis  method  uses  Equation  3.7.  The  damping  constant  and 
peak  pressure  (temporal)  method  uses  Equation  3. 10. 
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A.  ENERGY  using  the  THEORETICAL  DAMPING 
CONSTANT  and  PEAK  PRESSURE 

The  energy  contained  in  dominant  (primary)  bubbles,  as  defined  in  earlier  work 
[Jacobus,  1991;  Ostwald,  1992],  was  calculated  for  the  smooth  and  roughened  surface  data 
by  using  Equation  3.10.  Previously,  the  energy  of  only  the  dominant  bubble  of  each 
series  was  calculated  because  Equation  3.4  requires  knowledge  of  the  angle  0  and  a  range 
r  from  the  hydrophone.  Constrained  by  the  initial  definition  of  "primary  bubble",  only  the 
range  of  the  initial  Type  II  bubble  is  known.  Multiple  hydrophones  will  be  needed  to 
accurately  measure  the  range  to  all  bubbles  created.  However,  the  use  of  the  damping 
constant  and  the  peak  pressure  to  evaluate  the  energy  for  each  bubble  allows  for  an 
estimate  of  the  energy  in  the  entire  series. 

1.  Energy  at  Higher  Frequencies 

Figure  4. 1  (top)  shows  the  same  general  distribution  as  reported  by  Ostwald. 
Figure  4.1  (bottom)  shows  that  for  a  sloped  surface  at  the  time  of  impact,  there  exists 
more  energy  at  higher  frequencies.  Higher  frequency  bubbles  are  produced  by  smaller 
raindrops  [Kurgan,  1989;  Jacobus,  1991;  Ostwald,  1992].  Because  they  are  higher 
frequency,  one  would  assume  that  they  are  formed  by  the  Type  I  mechanism  (or  late 
aerosol  impacts).  There  is  no  definitive  proof  that  this  is  the  case.  The  roughened 
surface  could  possibly  produce  Type  II  bubbles  with  frequencies  higher  than  previously 
observed. 

2.  Time  Gap 

As  previously  stated,  earlier  work  used  a  shorter  time  series.  When  the 
duration  of  the  series  is  extended,  a  noticeable  time  gap  was  observed  for  the  smooth 
surface  data.  This  gap  occurs  between  1 10  and  140  ms  as  suggested  in  Figure  4.2  (top). 
This  time  gap  is  not  present  in  the  roughened  surface  data  [Figure  4.2,  bottom]. 
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Figure  4.1.  Comparison  of  Energies  of  Dominant  Bubbles  for  Smooth  (top)  and  Roughened 
(bottom)  sur&ces. 
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Energy  [Joules]  Energy  [Joules] 


(Dominont  Bubbles)  Energy  vs  Time  '^fter  Impoct.. .Smooth  Sfc 


(Oominont  Bubbles)  Energy  vs  Time  ofter  Impact. ..Roughened  Stc 


Figure  4.2.  Comparison  of  Dominant  Bubble  Energy  vs  Time  after  Impact  for 
Smooth  (top)  and  Roughened  (bottom)  surfaces. 


19 


When  the  frequencies  of  all  bubbles  that  were  present  in  the  time  series  (primary 
and  secondary)  are  plotted  against  the  time  after  impact,  the  time  gap  becomes  more 
prevalent  in  the  smooth  surface  data  [Figure  4.3,  top].  The  end  of  the  gap  (approximately 
140  ms  after  impact  of  the  large  drop)  corresponds  to  the  time  delay  for  an  aerosol 
trajectory  with  a  maximum  height  of  2.5  cm.  Bubbles  occurring  after  140  ms  would  not 
have  been  considered  as  dominant  bubbles  in  previous  studies.  Because  of  the  timing  as 
shown  in  Figure  3.1,  they  can  not  be  Type  II  bubbles. 

Aerosols  generated  by  impacting  a  smooth,  stationary  surface  tend  to  be 
"launched"  almost  vertically  and  have  been  seen  reaching  extreme  heights  of  40  cm  above 
the  surface.  This  corresponds  to  a  time  of  flight  (time  after  impact)  of  571  ms.  Visually, 
the  average  height  was  between  12  and  15  cm.  These  heights  correspond  to  times  after 
impact  of  312  and  350  ms,  respectively.  This  suggests  that  the  duration  of  the  data 
collection  should  be  lengthened  further,  although  there  would  be  a  corresponding 
reduction  in  sampling  rate. 

For  the  roughened  surface  data  in  Figures  4.2  and  4.3  (bottom)  the  time  gap  is 
not  present.  This  is  consistent  with  a  reduced  time  of  flight  caused  by  a  lower  trajectory 
of  the  aerosols.  The  rough  surface  aerosols  are  also  launched  with  a  much  greater 
horizontal  velocity  component  than  the  smooth  surface  aerosols.  It  was  quite  obvious  that 
the  horizontal  velocity  of  the  surface  wave,  as  well  as  the  wave  slope,  greatly  affect  the 
velocity  components  and  the  launch  angles  of  these  aerosols. 
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Figure  4.3  Comparison  of  Frequencies  of  ALL  Bubbles  vs  Time  after  Impact  for 
Smooth  (top)  and  Rough  (bottom)  sur&ces. 
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B.  SPECTRAL  ANALYSIS 

1.  Smooth  vs  Roughened  Surfaces 

Figure  4.4  compares  the  average  spectral  energy  density  of  the  smooth  surface 
data  (top  half)  to  the  average  spectral  energy  density  of  the  roughened  surface  data 
(bottom  half).  The  instantaneous  surface  slope  at  the  time  of  impact  was  not  measured 
because  we  were  looking  for  a  general  effect  due  to  a  randomly  sloped  surface.  It  can  be 
seen  from  Figure  4.4  that  the  average  total  energies  are  comparable.  This  was  not 
expected  because  more  high  frequency  bubbles  were  observed  in  the  roughened  surface 
data  compared  to  the  smooth  surface  data.  Typically,  high  frequency  bubbles  radiate  less 
energy  than  low  frequency  bubbles;  therefore  their  influence  on  the  total  energy  is  not  as 
significant  as  lower  frequency  bubbles. 

The  plot  of  the  average  energy  density  for  the  roughened  surface  in  Figure  4.4 
(bottom)  also  demonstrates  an  apparent  shift  of  the  peak  of  the  spectrum  to  higher 
frequencies.  This  was  expected  because  of  the  presence  of  higher  frequency  bubbles 
occurring  late  in  the  time  series. 
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Avq  Enerqy  Density  for  49  drops  for  Smooth  Sfc 
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Avq  Enerqy  Density  for  97  drops  for  Rouqhened  Sfc 
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Figure  4.4.  Comparison  of  Average  Energy  Densities  for  Smooth  (top)  and  Roughened 
(bottom)  surfaces. 


2.  Roughened  Surface  Categorized  by  Slope 

The  use  of  a  video  camera  allowed  us  to  determine  the  slope  of  the  surface  at  the 
time  of  impact.  The  slopes  were  grouped  into  three  categories.  The  mean  slopes  were  0, 
6,  and  11°  from  the  horizontal  (±  3°).  A  0°  slope  for  a  roughened  surface  is  different 
from  a  smooth  surface  because  of  the  horizontal  velocity  of  the  surface  waves.  The  wave 
velocity  was  seen  to  be  imparted  to  the  raindrop  at  impact  and  affected  the  nature  of  the 
aerosols  generated  by  the  impact. 

Figure  4.5  shows  the  three  slope  categories.  It  can  be  seen  that  the  average  total 
energies  are  again  comparable,  independent  of  the  slope  of  the  surface  at  the  time  of 
impact.  Forty-eight  drops  were  used  for  this  analysis.  The  energies  are  all  within  one 
standard  deviation  of  each  other.  Better  statistics  are  needed  before  stronger  statements 
can  be  made. 

As  with  the  randomly  roughened  surface  of  Figure  4.4,  the  peaks  of  the  spectra  of 
Figure  4.5  also  shift  to  higher  frequencies  as  the  slope  of  the  sur&ce  increases. 
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V.  CONCLUSIONS 


The  original  purpose  of  this  thesis  was  to  determine  the  effects  of  a  sloped  water 
surface  at  the  time  of  impact  for  a  4.6  mm  diameter  raindrop.  While  several  effects  were 
observed,  this  experiment  also  opens  the  door  to  numerous  other  studies.  The  results  of 
this  study  reveal  that; 

4  For  smooth  and  rough  surfaces  (including  slope  categories);  average  energies  per 
drop  are  of  the  same  order  of  magnitude  [Figure  4.4  and  Figure  4.5]. 

4  The  peak  of  the  energy  spectrum  for  the  sloped  surface  apparently  shifts  to  higher 
frequencies  [Figure  4.4,  bottom  and  Figure  4.5]. 

4  The  "time  gap"  discovered  for  smooth  surface  data  is  absent  for  the  sloped  surface 
data.  We  now  understand  that  the  time  gap  is  absent  in  the  sloped  surface  data 
because  the  characteristics  of  sloped  surface  aerosols  are  greatly  affected  by  the 
slope  of  the  surface  and  the  horizontal  velocity  of  the  surface  wave. 

Only  after  more  is  known  about  the  physics  of  underwater  sound  produced  by 
raindrops  of  different  sizes  generated  in  a  laboratory,  will  we  be  able  to  understand  the 
physics  of  underwater  sound  produced  by  natural  rainfall  at  sea.  This  will  ultimately 
establish  the  ability  to  estimate  rainfall  rate  by  means  of  remote  underwater  listening 
devices. 
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VI.  RECOMMENDATIONS 

A.  DEFINITIONS 

1.  Primary 

Redefine  Primary  (dominant)  to  be  the  most  energetic  bubble.  This  will 
remove  any  dependence  on  its  origin  (Type  I  or  H). 

2.  Type  I  Mechanism 

Only  include  bubbles  created  by  pinch-ofif  from  conical  cavity. 

3.  Type  II  Mechanism 

Only  include  bubbles  created  by  downward  jet  within  the  canopy  of  large  drop 

splashes. 

B.  FUTURE  EXPERIMENTS 

1.  More  Data  for  this  Experiment 

In  order  to  establish  respectable  statistics,  more  data  needs  to  be  taken.  More 
data  will  also  help  to  enforce  the  results  of  this  work.  The  effect  of  a  sloped  surface  for 
different  raindrop  sizes  should  be  studied. 

2.  Multiple  Raindrops 

So  far  only  single  raindrops  have  been  used  for  all  experiments.  Design  a 
multiple  drop  experiment  to  determine  if  there  is  any  interaction  between  the  drops  at 
impact  (linear  or  nonlinear).  The  number  and  volume  of  the  drops  must  be  controllable  in 
order  to  better  understand  the  resulting  spectra. 
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3.  Multiple  Hydrophones 

Multiple  hydrophones  can  be  used  to  determine  ranges  of  impacts  or  bubbles. 
This  will  determine  the  origin  (bubble  production  mechanism)  of  all  bubbles  occurring 
later  in  time. 

4.  Capillary  waves 

Generate  a  roughened  surface  by  creating  capillary  waves  and  then  perform 
this  same  experiment.  Compare  the  results  of  the  smooth,  sloped,  and  roughened 
(capillary)  surface  data. 

5.  Eliminate  the  %  Correction  Factor 

Record  the  data  series  even  if  bubbles  are  not  present.  This  will  allow  for  a 
more  accurate  average  total  energy  for  a  known  number  of  drops  and  bubbles. 

6.  Extend  the  Listening  Time 

Use  time  delays  of  320  ms  or  640  ms  to  extend  the  time  series  to  640  ms  or 
960  ms,  respectively.  This  will  not  affect  the  sampling  rate. 
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APPENDIX  A 


The  following  three  programs  were  written  to  analyze  the  single  drop  data  collected 
in  the  NFS  raindrop  tank.  The  data  was  collected  using  ComputerScope®  and  stored 
an  ASCn  file  format.  These  files  were  stripped  of  extraneous  headers  and  time 
information  and  stored  as  data  files  containing  only  the  output  voltage.  (The  time  step  is 
known  and  is  input  inside  these  programs).  The  stripped  ASCII  files  are  imported  into 
MATLAB.  These  programs  are  written  in  the  MATLAB  processing  language 
(MATLAB  *.m  files).  All  plotting  routines  can  be  easily  changed  within  the  program  to 
present  data  any  way  the  user  desires. 

The  first  program,  ANAL.M,  extracts: 

the  frequency  of  any  bubble  existing  in  the  time  series, 

the  time  after  impact  of  each  bubble, 

and  the  peak  to  peak  voltage  of  each  bubble  in  the  time  series. 

The  data  is  broken  up  into  "dominant”  and  "all".  It  then  saves  the  extracted  data  in 
matrices  used  in  succeeding  programs. 


%  Program  ANAL.M 

%  Program  to  ANALyze  data  taken  by  ComputerScope  and  converted  to  ASCII. 
%  Calculates  the  frequency,  time  after  impact,  and  Vpp  of  a  particular  bubble. 
%  Two  mouse  "ginputs”  are  required  by  the  user.  This  is  used  for  the  freq. 

%  It  also  creates  and  saves  (user  chooses  name  of  file)  the  data  in  a  *.mat  file 

%  for  further  manipulation  in  *.m  files  like  delta.m. 

clc 

clear 

clg 

dispfln  what  directory  will  I  find  the  data  ') 

direct  =  input('{i.e.  fiftl  path  [in  single  quotes]}  ? ....’); 

disp(") 

disp(") 

start  =  inputflnput  number  of  starting  file...'); 
stop  =  inputflnput  number  of  ending  file...'); 
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maxnum  =  stop  •  start  +  1; 

ftV=[  1;  %  empty  matrix  to  allow  the  building  of  the  true  data  matrix 

ftVdom={];  %  empty  matrix  to  allow  the  building  of  the  true  data  matrix 

clc 

disp('With  what  letter  does  the  data  file  name  start  ') 
a=input('{i.e.  r*.dat,  s*.dat,  etc}  [in  single  quotesj...?'); 
for  i=  limaxnum; 
clg 
clc 

eval(['load  ',direct,T,a,int2str(start  +  i  - 1),  '.dat']) 
eval(['v  =  ',a,int2str(start  +  i  - 1),';']) 
eval(['clear  ',a,int2str(start  +  i  -  l),';'j) 
delt=20e-6; 

V  =  V  -  mean(v); 
plot(v) 
xlabel('n') 
ylabel('[volts]') 

title(['Drop  #  ',num2str(i),'  for  this  data  set.']) 

grid 

pause 

clc 

disp('Entering  keyboard  mode  so  that  the  user  can  determine  the  #  of) 
dispCbubbles  for  this  time  series  (in  case  of  uncertainty).*) 
disp(") 

di^'[CTRL-zl  to  exit  k'^  br  .rd  mode') 

dispC') 

disp{") 

keyboard 

clc 

bub=rinput('How  many  bubbles  in  this  series?...'); 

forh=l:bub; 

clc 

low=input('Input  lower  data  point  (n)  for  plotting...'); 

up=input('Input  upper  data  point  (n)  for  plotting...'); 

subplot(211),plot(v) 

xlaWCn*) 

ylabcl('[volts]') 

title(['Drop  #  ',num2str(i),'  for  this  data  set']) 
grid 

sulq)Iot(2 1 2),plot(v(low;up)) 

xlatelCn*) 

ylabcl('[volts]') 

title(['Bubble  #  ',num2str(bub),'  for  this  drop  ]) 

grid 

pause 

disp(") 

check«input('Are  these  the  correct  limits?... 1-Yes  2-No...’); 
while  check  =  2 
clc 

disp(") 
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disp('Last  limits  were... ') 
dispC’) 

liinits={low  up] 
disp(") 

iow=input('Input  lower  data  point  for  plotting...'); 
up=nnput('Input  upper  data  point  for  plotting...'); 
clg 

subplot(211),plot(v) 

xlaM('n') 

ylabel('(volts]') 

title(['Drop  #  ',num2str(i),'  for  this  data  set']) 
grid 

subplot(2 12),plot(v(low;up)) 

xlatel('n') 

ylabel('[volts]') 

title(['Bubble  #  ',nuin2str(bub),'  for  this  drop']) 

grid 

pause 

clear  check 

check=input('Are  these  the  correct  limits?...  1-Yes  2-No...'); 
end  %  end  while 
clc 

[mx,c]=max(v(low;up)); 

mn=min(v(low:up)); 

[xtyv]=ginput(2); 

freq=(abs(l/((xt(l).xt(2))*delt)))/1000  %  [kHz] 
t*((c+low)*delt)*1000  %[ms] 

Vpp=abs(mx)  +  abs(mn)  %  [volts] 

disp(") 

flVb=[freq  t  Vpp];  %  freq,  time,  Vpp  for  each  buU>le 

disp(") 

ifbub~=l; 

dom=input('Is  this  the  dominant  bubble  of  the  series?...  1-Yes  2-No...'); 
else 
dom=l; 
end 

if  dom==l; 

ftVdom=(ftVdom;ftVb]; 

end 

ftV=[ftV;ftVb]; 

clg 

plot(v) 

title(['Extra  look  at  drop  #  ',num2str(i),'  for  this  data  set']) 

grid 

pause 

clc 

end  %  next  b 
clc 

clear  V 

end  %  next  i 

mat=input('What  name  for  the  *.mat  file  [in  single  quotes]  ? ....'); 
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dispC'  ’) 

disp(The  matrix  containing... freq  [kHz]  t  [ms]  Vpp  [volts)...«ill  be  saved') 
disp('for  future  uses  (i.e.  delta.m).  The  variables  will  be  ftV  and  ftVdom.') 
disp(") 

disp('lnsert  a  disk  into  drive  a:  to  receive  the  ♦.mat  file...') 

disp('[ENTER]  to  continue') 

pause 

!a: 

eval(['save  '.mat,'  ftV  flVdom']) 

!c: 

clc 

disp('It  is  now  recommended  to  set  up  and  nm  DELTA.M') 
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The  second  program,  DELTA.M,  uses  the  matrices  formed  in  ANAL.M  and 
calculates  the  theoretical  damping  constant  (5).  Then  the  energy  using  5  and  the  1/2  of 
Vpp  ( peak  voltage)  found  in  ANAL.M  is  calculated.  Finally,  many  different  combinations 
of  data  are  plotted  like  Figures  4. 1,  4.2  and  4.3. 


%  Program  DELTA.M 

%  to  be  used  after  running  ANAL.M  and  before  ANAL2.M 

%  Creates  monster  matrix  and  calculates  the  damping  constants  and 
%  Energy  based  on  these  deltas  (Clay  &  Medwin  text  Ch  6) 

%  Plots  different  combinations  of  the  data. 

format  short  e 

clc 

clear 

pack 

clg 

dispCData  consists  of  the  matrices  created  in  ANAL.M.  It  is’) 
dispCthe  *.mat  file  containing  the  matrices  consisting  of;*) 
dispC ') 

dispCl  f  (kHz)  t(ms)  Vpp  (volts)  ].') 

dispC  ’) 

disp(") 

dispC  What  was  the  name  of  the  '*.mat  file  from') 

file  =  inputCANAL.M  {FULL  PATH  including  filename  (in  single  quotes)} . ?'); 

clc 

dispCWhere  are  the  files  d.dat  and  z.dal') 

direct2  =  input('{i.e.  FULL  PATH  (in  single  quotes)} . ?'); 

eval(['load  ',file]); 
evaI(['load  ',direct2,y,’d.dat']); 
eval(['load  ',direct2,T,'z.dat']); 


d=d./100;  %  [m] 

z=z./100;  %  [m] 

d=d(l:length(ftVdom(:,l)));  %[m] 

z=z(l:length(ftVdom(:,l)));  %[m] 


%  rhoc  =  inputflnput  the  Specific  Acoustic  Impedance  (rhoc)...'); 
%  diam  =  inputflnput  drop  diameter  [mm]...'); 

%  gain  =  inpuU'Input  hydrophone  gain...'); 

%  ML  =  inputflnput  hydrophone  Sensitivity  [V/Pa]...'); 
rO=l;  %  1  meter  on  axis 
c  =  1500;  %  sound  speed 

rhoc  =  L54e06;  %  for  seawater 

ML  ■=  2.66C-5;  %  Hydrophone  sensitivity  (V/Pa)  over  1-20  kHz 
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diam=4.6; 

gain=100; 


%  calculate  the  damping  constant  based  on  a  linear  regression  of  the 
%  delta  vs.  freq  (log-log  plot)  from  Clay  &  Medwin 

del=.025  ♦  ftVdom(:, I).''!  1/3);  %  ftVdom(:,  1)  is  freq  column  [kHz] 

%  C^culate  energy  based  on  damping  constant  'delta* 

r=sqrt{z.''2  +  d.''2);  %  [m] 

k=2'''pi*ftVdom(:,  1)*  1000/c;  %  1000  to  convert  to  Hz 

CF=(r.^2)./(gain*ML*rO*z); 

NFC  =  (k.*r) ./  (sqrt((k.*r).^  2  +  1)); 

pax=NFC.*CF.'*ftVdom(:,3)/2;  %  1  m  on  axis  based  on  PEAK  voltage 
E=((r.*pax).''‘2)./(1000*6*rhoc.*ftVdom(;,l).’*del);  %  1000  to  convert  to  Hz 

ftVdEdom=[ftVdom  del  E  ]; 

%  This  part  of  delta.m  plots  different  combinations  of  ftVdEdom. 

%  This  matrix  now  consists  of  five  columns; 

%  [f(kHz)  t(msec)  Vpp(Volts)  delta(damping  constant)  Energy(joules)] 
clc 

dispCWhat  type  of  surface  was  used  for  this  data') 

surf  =  input('{i.e.  Smooth  or  Roughened  (in  single  quotes)} . ?'); 

clc 

name  =  input('What  name  do  you  want  for  the  •.met  file  {in  single  quotes} . ?  ’) 

clc 

disp('What  name  do  you  want  for  the  •.mat  file  containing ') 

name2  =  input('the  ( f  t  Vpp  delta  E  ]  matrix  {in  single  quotes} . ?  *); 

clc 

m=mean(flVdEdom); 

%plot#l 

%  Dominant  bubble  Energy  vs  freq 
%  11={0,  25, -12, -6]; 

%  axis(ll) 

%  subplot(211),semilogy(ftVdEdom(:,l),flVdEdom(:,5),'+') 

%  tit=['Dominant  Energy  vs  Frequency.. .',surf,'  Sfc']; 

%  title(tit) 

%  xlabeK'freq  [kHz]') 

%  ylabelCEnergy  [Joules]') 

%  gtext(['Avg  Energy  =  ',num2str(m(5)),'  Joules']) 

%  plot  #2 

%  Dominant  bubble  Energy  vs  time 
%  12=[0,  300, -12, -6]; 

%  axis(12); 

%  subplot(2 1  l),semilogy(flVdEdom(:,2),flVdEdom(;,5),'+') 

%  sul^lot(2 12),semilogy(ftVdEdom(:,2),flVdEdom(:,5),'+') 

%  tit^CDominant  Energy  vs  Time  after  Impact. .  .',surf,'  Sfc'] ; 
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%  title(tit) 

%  xlabelCtime  [msec]') 

%  ylabelCEnergy  [Joules]’) 

%clc 
%  pause 

%  disp('Insert  a  disk  into  drive  a  to  receive  the  ♦.met  and  *.mat  files.') 

%  dispCENTER  to  coninue') 

%  pause 
%  !a: 

%  eval(('meta  ',namel)  %  plots  #  1  &  #2 
%  !c: 

%  clg 

%  plot  #3 

%  Dominant  bubble  Energy  vs  drop  # 

%  subplot(2 1  l),semilogy(ftVdEdom(:,5),'+') 

%  tit=('Dominant  Energy  vs  Drop  number...’, surf,’  Sfc’]; 

%  title(tit) 

%  xlabeU’Drop  #’) 

%  ylabeU’Energy  [Joules]’) 

%  plot  #4 

%  All  bubbles... Freq  vs  time 

14=[0,  300, 0,  25];  %  force  the  axis  limits  for  direct  comparisons 

axis(14) 

subplot(21  l),plot(ftV(:,2),ftV(;,  1),’+’) 
gtextCALL  bubbles') 

%  subplot(212),plot(ftV(;,2),ftV(:,l),’+’) 

tit=[’Frequency  vs  Time  after  Impact...  ’.surf,’  Sfc’]; 
title(tit) 

xlabeK’t  [msec]’) 
ylabeU’freq  [kHz]’) 

pause 
%  !a: 

eval([’meta  ’.name]) 
eval([’save  ’.name,’  ftVdEdom’]) 

%  !c; 
clc 

disp(’  ’) 
disp(’  ’) 

disp(’It  is  now  recommended  that  you  set  up  and  run  ANAL2.M') 
disp(’  ’) 

dispCThank  you  for  playing  DELTA.’) 
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The  final  program,  ANAL2.M,  calculates  and  averages  spectra  for  a  series  of  drops 
and  produces  plots  like  those  shown  in  Figures  4.4  and  4.5. 

%  Program  ANAL2.M 
%  for  rough  and  smooth  surfaces 

%  program  to  ANALyze  spectral  level  determination  of  raindrop  data 
%  Converts  bubble  time  series  to  power  spectrum.  Outputs  result. 

clc 
clear 
clg 

c  =  1500; 
rhoc  =  1.54e06; 

ML  =  2.66e-5; 
diam=4.6; 
gain=100; 

N=16384; 
delt  =  20  •  10''(-6); 
dispC  ’) 

dispC  NOTE:  Drop  file  names  must  be  of  the  form  r#.dat  or  s#.daL') 
dispC  The  drop  files  must  be  numbered  consecutively.') 
disp(") 

dispC  Ensure  impaa  distances  are  in  file  d.dat. ') 
dispC  and  hydrophone  depths  are  in  file  z.dat.') 

dispC  D.DAT  and  Z.DAT  MUST  be  in  the  same  directory  as  your  data  files.’) 
dispC  Distances  must  be  in  a  single  column,  in  cm') 
dispC ') 


%  sound  speed  (use  the  conect  sound  speed  for  the  salinity) 
%  for  seawater 

%  Hydrophone  sensitivity  (V/Pa)  over  1-20  kHz 


dcheck  =  inputCDid  you  update  d.dat  (distance)  and  z.dat  (depth)  files  (1  =  yes,  2  =  no)..’); 
if  dcheck  1 
dispC') 

errorCUpdate  the  files  and  rerun  this  program.') 
end 

clear  dcheck; 
clc 

%  diam  =  input('Input  drop  diameter  (mm)..,'); 

%  gain  =  input('Input  hydrophone  gain...'); 

%  ML  =  inputCInput  hydrophone  Sensitivity  [V/Pa]...'); 

%  N  =  inputCInput  number  of  points  for  FFT...'); 

%  delt  =  inputCInput  sample  period  (time  between  samples,  usee)...'); 

%  delt  =  delt  •  10^(-6); 

pack 

clc 
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surf  =  inputC  What  type  of  surface?  (1 -Smooth  l-Roughened)  '); 
dispC  ’) 


if  surf=2  %  Rough  Surfaces 
clc 

disp('DATA  includes  z.dat,  d.dat.  and  r^.dat  files') 
disp(”) 

disp('In  what  directory  will  I  find  the  DATA') 

direct  =  input('{i.e.  fiill  path  [in  single  quotes]}  ? ...'); 

eval(['load  ',dir^ , '\d.dat;']) 

d=d./100;  %[m] 

eval(['load  '.direct  .^.dat;']) 

z=z./100;  %[m] 

clc 

name  =  input('What  do  you  want  to  call  the  meta  file  [in  single  quotes]  ? ...’); 
disp(") 

disp('What  is  the  Sfc  slope  for  this  data  (i.e.  "5.5  degrees")') 

slope  =  input('[in  single  quotes]  ?...'); 

clc 

dispCFor  Roughened  Surfaces...') 
start  =  inputflnput  number  of  starting  file...'); 
stop  =  inputflnput  number  of  ending  file...'); 
maxnum  =  stop  -  start  +  1; 

perc  =  inputCInput  %  of  drops  that  produce  at  least  one  bubble...'); 
perc  =  perc  / 100; 
df  =  1/(N  •  delt); 

Esum  =zeros(  1  :(N/2));  %  Initializing  size  of  row  vector 

(df.*  [0:l:(N/2)-l])'; 
k  =  (2*pi.*f)./c; 

iO=l;  %  [mj 

for  i  =  l:maxnum; 

eval(['Ioad  ',direct,'\r*,int2str(start  +  i  -  1),  '.dat']) 
eval(['v  =  ','r',int2str(start  +  i  -  1),';']) 
eval(('clear  ',’r',int2su(start  +  i  - 1),';']) 

V  =  V  -  mean(v); 

%  t=0:delt:(length(v)-l)*delt; 

V  =  abs(ffl(v,N)); 
clear  V 

r(i)  =  sqrt(z(i)''2+d(i).''2); 

NFC  =  (k-^rfi)) ./  (sqrt((k.'»r(i)).''  2  +  1)); 
CF=((perc*r(i)^2)/(gain*ML*z(i)*iO)); 

Pax  =  CF  .*  NFC  .•  V(l:(N/2)); 

%  Using  eqn  41  from  thesis  by  Ostwald 
EH(2*delt)./(N*d0).*  Pax(l:N/2).''2; 
clear  Pax  V 

E=((2'»pi*r(i).''2)/(3*rhoc))*Ef; 
clear  Ef 

Esum  =  Esum  +  E'; 
clear  E 
end 


%  [m] 

%  [Pa/v] 

%  [Pa]  1  m  on  axis 
%  [Pa''2-s/Hz] 

%  [Joules/Hz] 

%  for  averaging  purposes  later 
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clear  CF  NFC  start  stop  N  delt  i  surf  r 


Eavg  =  Esum/maxnum;  %  Average  Energy  density  of  all  drops 

etot=sum(Eavg)*df;  %  Evaluate  the  integral  eqn  42 

SL  =  (10*logl0(Esum)  +  120)';  %  dB  re  1  uPa'^2/Hz 
clc 

dispCFour  gtexts  follow;  Sfc  slope,  Avg  Energy,  slope  dependence,  freq  shift') 

disp('[Enter  to  continue]') 

pause 

I1=(0, 18, 0, 4e-13]; 
axisfll) 

%  subplot(2 1  l),plot(f./1000,Eavg); 

%  subplot(212),plot(f./1000,Eavg); 
plot(t/1000,Eavg); 
xlabeK'freq  (kHz)') 
ylabel('[Joules/Hz]') 
if  maxnum=l; 

title(['Avg  Energy  Density  for  ',num2str(niaxnum),'  drop  for  Roughened  Sfc']) 
else 

titled'Avg  Energy  Density  for  ',num2str(maxnum),'  drops  for  Roughened  Sfc']) 
end 

%({'Drop  Diameter  =  ',num2str(diam),'  mm']) 
gtext(('Sfc  slope  =  ',slope]) 
gtextff  Average  Total  Energy  =  ',num2str(etot),'  [J]']) 
gtextCNo  significant  dependence  on  slope') 
gtextCSpectrum  shifts  to  higher  frequencies') 

%  subplot(212),plot(f(2:8190)./1000,SL(2:8l90)); 

%  xlateK'freq  (kHz)') 

%  ylabelCdBreluPa''2/Hz') 

%  ifmaxnum=l; 

%  title(['Spectrum  Level  for  ',num2str(maxnum),' drop']) 

%  else 

%  title(['Spectrum  Level  for  ',num2str(maxnum),' drops']) 

%  end 


pause 
%  clc 

%  disp('Insert  the  *.met  disk  into  A') 
%  disp('[Enter  to  continue]') 

%  pause 
%  !a; 

eval(['meta  ',name]) 

%  !c: 
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else  %  Smooth  Surfaces 


clc 

dispCDATA  includes  d.dat,  z.dat,  and  s*.dat  files') 

direct  =  input('In  what  directory  will  I  find  the  data  {i.e.  full  path  [in  single  quotes]}  ?  .... 

eval(['load  ',direct,'\d.dat']) 

eva](['load  ',direct,'\z.dat']) 

d=d./100;  %[m] 

z=z./100;  %  [m] 

name  =  input(’What  do  you  want  to  call  the  meta  file  [in  single  quotes]}  ? ...'); 
clc 

disp('For  Smooth  Surfaces...') 
start  =  inputflnput  number  of  starting  file...'); 
stop  =  inputflnput  number  of  ending  file...'); 
maxnum  =  stop  -  start  +  1; 

perc  =  inputf'Input  %  of  drops  that  produce  at  least  one  bubble...'); 
perc  =  perc  / 100; 
df  =  1/(N  •  delt); 

Esum  =zeros(  1  :G^/2));  %  Initializing  size  of  row  vector 

f=(df.*[0.1:(N/2)-l])'; 
k  =  (2*pi.'*f)./c; 

r  =  sqrtfz.^  2  +  d.''  2);  %  [m] 

rO=l;  %  [m] 

for  i  =  Irmaxnum; 

evaI(['load ', direct, '\s',int2str(start  +  i  - 1),  '.dat')) 
eval(['v  =  ','s',int2str(start  +  i  - !),';']) 
cval(('clear  ','s',int2str(start  +  i  - 1),';']) 

V  =  V  -  meanfv); 

%  t=0:delt;(length(v)-l)'''delt; 

V  =  abs(fift(v,N)); 
clear  V 

i  %  Progress  Pointer  for  the  user 
NFC  =  (k.'*r(i)) ./  (sqrt((k.*r(i)).''  2  +  1)); 

CF={(perc*r(i)^2)/(gain*ML*z(i)*rO)); 

Pax  =  CF  .*  NFC  .*  V(l:(N/2)); 

%  Using  eqn  4 1  from  thesis  by  Ostwald 
Ef=((2*dclt)./(N*df)).*  Pax(l:N/2).''2; 
clear  Pax  V 

E=((2*pi*r(i).''2)/(3*rhoc))*Ef; 
clear  Ef 

Esum  =  Esum  +  E'; 
clear  E 
end 

clear  CF  NFC  start  stop  N  delt  i  surf  r 

Eavg  ^  Esum/maxnum; 
etot*sum(Eavg)*df; 

SL  -  (10*logl0(Esum)  +  120)'; 


%  [Pa/v] 

%  [Pa]  1  m  on  axis 
%  [Pa''2-s/Hz] 

%  [Joules/Hz] 

%  for  averaging  purposes  later 


%  Average  Energy  density  of  all  drops 
%  Evaluate  the  integral  eqn  42 
%dBieluPa^2/Hz 
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disp(") 

%  disp('Please  insert  the  *.met  disk  into  drive  A  (Enter  to  continue]') 

%  pause 
clc 

dispCOne  gtext  follows:  Avg  Energy*) 

di^'[Enter  to  continue]') 

pause 

11=(0, 18,0,4e-13]; 
axis(ll) 

%  subplot(211),plot(f./1000,Eavg); 
plot(f/1000,Eavg); 
xiabel('freq  (kHz)') 
ylabel('[Joules  /  Hz]') 
if  inaxnuin=l; 

titlefCAvg  Energy  Density  for  '.numlstifmaxnuin),'  dix^  for  Smooth  Sfc']) 
else 

title(['Avg  Energy  Density  for  ',num2str(maxnum),'  drops  for  Smooth  Sfc']) 
end 

%gtext(['Dtop  Diameter=  ',num2str(diam),'  mm']) 
gtextCf  Average  Total  Energy=  ',num2str(etot),'  [J]*]) 

%  subplot(212),plot(fl[2;8190)./1000,SL(2:8190)); 

%  xlatel('freq  (kHz)') 

%  ylabel('dBreluPa''2/Hz') 

%  ifmaxnum— 1; 

%  title(['Spectrum  Level  for  ',num2str(maxnum),' drop']) 

%  else 

%  title(['Spectrum  Level  for  ',num2str(maxnum),’ drops’]) 

%  end 


pause 
%  !a: 

eval(['meta  '.name;]) 
%  !c: 
end 
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